---
title: "CARRS RDD analysis using a standardized cutoff"
knit: (function(input, ...) {rmarkdown::render(input, output_file="")})
output: pdf_document
geometry: "left=2cm,right=2cm,top=1cm,bottom=1cm"
date: "`r format(Sys.time(), '%d %B, %Y')`"
---

```{r , setup, include=FALSE, eval = FALSE}
knitr::opts_chunk$set(
  comment = '', warning=FALSE,  message=FALSE, results= FALSE , fig.width = 9,fig.height = 7, echo=FALSE, root.dir =""
  
)

knitr::opts_knit$set(
  root.dir =""
)

```

```{r packages}
rm(list=ls())
  library(tidyverse)
  library(gridExtra)
  library(dplyr)
  library(rdrobust)
  library(rddensity)
  library(viridisLite)
  library(foreign)
  library(ggplot2)
  library(data.table)
  library(doBy)
  library(socviz)
  library(openxlsx) 
  library(modelsummary)
  library(labelled)
  library(gtsummary)
  library(gt)
  library(kableExtra)
  library(ggpubr)
```

```{r loaddata}
data_treatment<-data.frame(read.csv("carrs_full.csv"))
data_diagnosis<-data.frame(read.csv("carrs_selected.csv"))
data_bp<-data.frame(read.csv("carrs_bp.csv"))
```

```{r categorical}
##Defining the common function##
  ctrls = function(a,b,d) {
    
    d$Y<-a
    d$X<-b

    out<- rdrobust(d$Y, d$X, c = 0, 
                    cluster =d$pid, 
                    covs = cbind(d$w01, d$w23, d$w45), 
                    kernel = "triangular", 
                    p = 1,rho=1)

    ret <- data.frame(
    sex = c(formatC(round(100*out$coef[1,1],2),2,format="f"), formatC(round(out$pv[3,1],2),2,format="f"), 
            paste0( "(", formatC(round(100*out$ci[3,1],2),2,format="f"), " - ", formatC(round(100*out$ci[3,2],2),2,format="f"), ")" ), 
            formatC(round(out$bws[1,1],2),2,format="f"), formatC(out$N_h[1] + out$N_h[2])))
    
    rownames(ret) <- c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")
  
    output<<-ret %>% add_rownames 
    out<<-out
  }
```

```{r continuous}
##Defining the common function##
  ctrlscont = function(a,b,d) {
    
    d$Y<-a
    d$X<-b

    out<- rdrobust(d$Y, d$X, c = 0, 
                    cluster =d$pid, 
                    covs = cbind(d$w01, d$w23, d$w45), 
                    kernel = "triangular", 
                    p = 1,rho=1)

    ret <- data.frame(
    sex = c(formatC(round(out$coef[1,1],2),2,format="f"), formatC(round(out$pv[3,1],2),2,format="f"), 
            paste0( "(", formatC(round(out$ci[3,1],2),2,format="f"), " - ", formatC(round(out$ci[3,2],2),2,format="f"), ")" ), 
            formatC(round(out$bws[1,1],2),2,format="f"), formatC(out$N_h[1] + out$N_h[2])))
    
    rownames(ret) <- c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")
  
    output<<-ret %>% add_rownames 
    out<<-out
  }
```

**Placebo outcome: Sex**

```{r female}
data_diagnosis<-as.data.table(data_diagnosis)
data_treatment<-as.data.table(data_treatment)
data_bp<-as.data.table(data_bp)

ctrls(data_diagnosis$female,data_diagnosis$runningvar,data_diagnosis)
rob_femaleD <-output%>%mutate(Total=sex)%>%select(-sex)
summary(out)

ctrls(data_treatment$female,data_treatment$runningvar,data_treatment)
rob_femaleT <-output%>%mutate(Total=sex)%>%select(-sex)
summary(out)

ctrls(data_bp$female,data_bp$runningvar,data_bp)
rob_femaleB <-output%>%mutate(Total=sex)%>%select(-sex)
summary(out)

```

**Placebo outcome: Education**

```{r edu}

ctrls(data_diagnosis$edu,data_diagnosis$runningvar,data_diagnosis)
rob_eduD <-output%>%mutate(Total=sex)%>%select(-sex)
summary(out)

ctrls(data_treatment$edu,data_treatment$runningvar,data_treatment)
rob_eduT <-output%>%mutate(Total=sex)%>%select(-sex)
summary(out)

ctrls(data_bp$edu,data_bp$runningvar,data_bp)
rob_eduB <-output%>%mutate(Total=sex)%>%select(-sex)
summary(out)

```

**Placebo outcome: Age**

```{r age}

ctrlscont(data_diagnosis$age_calc,data_diagnosis$runningvar,data_diagnosis)
rob_ageD <-output%>%mutate(Total=sex)%>%select(-sex)
summary(out)

ctrlscont(data_treatment$age_calc,data_treatment$runningvar,data_treatment)
rob_ageT <-output%>%mutate(Total=sex)%>%select(-sex)
summary(out)

ctrlscont(data_bp$age_calc,data_bp$runningvar,data_bp)
rob_ageB <-output%>%mutate(Total=sex)%>%select(-sex)
summary(out)
```

**Placebo outcome: Religion**

```{r religion}
data_diagnosis <- data_diagnosis %>% 
  mutate(religion_num = as.numeric(factor(religion)))
data_treatment <- data_treatment %>% 
  mutate(religion_num = as.numeric(factor(religion)))
data_bp <- data_bp %>% 
  mutate(religion_num = as.numeric(factor(religion)))

ctrls(data_diagnosis$religion_num,data_diagnosis$runningvar,data_diagnosis)
rob_religionD <-output%>%mutate(Total=sex)%>%select(-sex)
summary(out)

ctrls(data_treatment$religion_num,data_treatment$runningvar,data_treatment)
rob_religionT <-output%>%mutate(Total=sex)%>%select(-sex)
summary(out)

ctrls(data_bp$religion_num,data_bp$runningvar,data_bp)
rob_religionB <-output%>%mutate(Total=sex)%>%select(-sex)
summary(out)

```

**Placebo outcome: Diabetes diagnosis**

```{r diabetes}

ctrls(data_diagnosis$diabetesknow,data_diagnosis$runningvar,data_diagnosis)
rob_diabD <-output%>%mutate(Total=sex)%>%select(-sex)
summary(out)

ctrls(data_treatment$diabetesknow,data_treatment$runningvar,data_treatment)
rob_diabT <-output%>%mutate(Total=sex)%>%select(-sex)
summary(out)

ctrls(data_bp$diabetesknow,data_bp$runningvar,data_bp)
rob_diabB <-output%>%mutate(Total=sex)%>%select(-sex)
summary(out)

```

**Placebo outcome: Heart attack**

```{r heart}

ctrls(data_diagnosis$heartknow,data_diagnosis$runningvar,data_diagnosis)
rob_heartD <-output%>%mutate(Total=sex)%>%select(-sex)
summary(out)

ctrls(data_treatment$heartknow,data_treatment$runningvar,data_treatment)
rob_heartT <-output%>%mutate(Total=sex)%>%select(-sex)
summary(out)

ctrls(data_bp$heartknow,data_bp$runningvar,data_bp)
rob_heartB <-output%>%mutate(Total=sex)%>%select(-sex)
summary(out)

```

**Placebo outcome: Stroke**

```{r stroke}

ctrls(data_diagnosis$strokeknow,data_diagnosis$runningvar,data_diagnosis)
rob_strokeD <-output%>%mutate(Total=sex)%>%select(-sex)
summary(out)

ctrls(data_treatment$strokeknow,data_treatment$runningvar,data_treatment)
rob_strokeT <-output%>%mutate(Total=sex)%>%select(-sex)
summary(out)

ctrls(data_bp$strokeknow,data_bp$runningvar,data_bp)
rob_strokeB <-output%>%mutate(Total=sex)%>%select(-sex)
summary(out)

```

**Tables**

```{r}
female_table<-merge(rob_femaleD,rob_femaleT, by="rowname")
female_table<-merge(female_table,rob_femaleB, by="rowname")

female_table<-female_table %>% arrange(match(rowname, c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")))

colnames(female_table)[2] <- "Diagnosis"   
colnames(female_table)[3] <- "Treatment"   
colnames(female_table)[4] <- "Blood pressure"   

age_table<-merge(rob_ageD,rob_ageT, by="rowname")
age_table<-merge(age_table,rob_ageB, by="rowname")
age_table<-age_table %>% arrange(match(rowname, c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")))

colnames(age_table)[2] <- "Diagnosis"   
colnames(age_table)[3] <- "Treatment"   
colnames(age_table)[4] <- "Blood pressure"   

religion_table<-merge(rob_religionD,rob_religionT, by="rowname")
religion_table<-merge(religion_table,rob_religionB, by="rowname")
religion_table<-religion_table %>% arrange(match(rowname, c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")))

colnames(religion_table)[2] <- "Diagnosis"   
colnames(religion_table)[3] <- "Treatment"   
colnames(religion_table)[4] <- "Blood pressure"   

educ_table<-merge(rob_eduD,rob_eduT, by="rowname")
educ_table<-merge(educ_table,rob_eduB, by="rowname")
educ_table<-educ_table %>% arrange(match(rowname, c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")))

colnames(educ_table)[2] <- "Diagnosis"   
colnames(educ_table)[3] <- "Treatment"   
colnames(educ_table)[4] <- "Blood pressure"   

diab_table<-merge(rob_diabD,rob_diabT, by="rowname")
diab_table<-merge(diab_table,rob_diabB, by="rowname")
diab_table<-diab_table %>% arrange(match(rowname, c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")))

colnames(diab_table)[2] <- "Diagnosis"   
colnames(diab_table)[3] <- "Treatment"   
colnames(diab_table)[4] <- "Blood pressure"   

heart_table<-merge(rob_heartD,rob_heartT, by="rowname")
heart_table<-merge(heart_table,rob_heartB, by="rowname")
heart_table<-heart_table %>% arrange(match(rowname, c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")))

colnames(heart_table)[2] <- "Diagnosis"   
colnames(heart_table)[3] <- "Treatment"   
colnames(heart_table)[4] <- "Blood pressure"   

stroke_table<-merge(rob_strokeD,rob_strokeT, by="rowname")
stroke_table<-merge(stroke_table,rob_strokeB, by="rowname")
stroke_table<-stroke_table %>% arrange(match(rowname, c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")))

colnames(stroke_table)[2] <- "Diagnosis"   
colnames(stroke_table)[3] <- "Treatment"   
colnames(stroke_table)[4] <- "Blood pressure"   

```
##table making function


```{r}
#combining diagnosis and treatment
  filename <- paste0("rob_ctrls.tex")
  caption=paste0("Pre-baseline control variables")

  comb_table<-rbind(female_table,age_table,religion_table,educ_table,diab_table,heart_table,stroke_table)

combined<-kbl(comb_table,caption=caption, booktabs = T, linesep = '', escape = FALSE,align = "lccccc", format = "latex",col.names = c(" ","Diagnosis","Treatment","Blood pressure"),digits = 2) %>%
    kable_styling(latex_options = "HOLD_position", position = "center") %>%
    pack_rows("Panel A: Female", 1,5) %>%
    pack_rows("Estimation results", 1,3) %>%
    pack_rows("Bandwidth details", 4, 5) %>%
    pack_rows(" ", 6,10) %>%
    pack_rows("Panel B: Age", 6,10) %>%
    pack_rows("Estimation results", 6, 8) %>%
    pack_rows("Bandwidth details", 9, 10) %>%
    pack_rows(" ", 11,15) %>%
    pack_rows("Panel C: Religion", 11,15) %>%
    pack_rows("Estimation results", 11, 13) %>%
    pack_rows("Bandwidth details", 14, 15) %>%
    pack_rows(" ", 16,20) %>%
    pack_rows("Panel D: Education", 16,20) %>%
    pack_rows("Estimation results", 16, 18) %>%
    pack_rows("Bandwidth details", 19, 20) %>%
    pack_rows(" ", 21,25) %>%
    pack_rows("Panel E: Diabetes", 21,25) %>%
    pack_rows("Estimation results", 21, 13) %>%
    pack_rows("Bandwidth details", 24, 25) %>%
    pack_rows(" ", 26,30) %>%
    pack_rows("Panel F: Heart attack", 26,30) %>%
    pack_rows("Estimation results", 26, 28) %>%
    pack_rows("Bandwidth details", 29, 30) %>%
    pack_rows(" ", 31,35) %>%
    pack_rows("Panel G: Stroke", 31,35) %>%
    pack_rows("Estimation results", 31, 33) %>%
    pack_rows("Bandwidth details", 34, 35) %>%
    footnote(general = c("xx","xx"),general_title="",threeparttable = T, fixed_small_size = T)%>%
    save_kable(filename, format = "latex")

combined
```